Hard ellipsoids: analytically approaching the exact overlap distance 
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Following previous work (JCP 134, 201103 (2011)), the replica exchange Monte Carlo technique 
is used to produce the equation of state of hard 1:5 aspect-ratio oblate ellipsoids for a wide density 
range. Here, in addition to the analytical approximation of the overlap distance given by Berne 
and Pechukas (BP) and the exact numerical solution of Perram and Wertheim, we tested a simple 
modification of the original BP approximation (MBP) which corrects the known T-shape mismatch 
of BP for all aspect-ratios. We found that the MBP equation of state shows a very good quantitative 
agreement with the exact solution. The MBP analytical expression allowed us to study size effects 
on the previously reported results. For the thermodynamic limit, we estimated the exact 1:5 hard 
ellipsoid isotropic-nematic transition at the volume fraction 0.343 ± 0.003, and the nematic-solid 
transition in the volume fraction interval (0.592 ± 0.006) - (0.634 ± 0.008). 

PACS numbers: 64.30.-t, 64.70.mf, 61.30.Cz 



I. INTRODUCTION 

In 1972 Berne and Pechukas [1] (BP) published a pio- 
neer work where molecules were represented by an uni- 
axial ellipsoidal Gaussian and their pair interaction was 
then related to their overlap integral. This was proba- 
bly one of the first attempts to produce a coarse-grained 
model and it certainly was a keystone for the develop- 
ment of the widely used anisotropic Gay-Berne [2[ inter- 
action. This last pair-potential is very popular since it 
allows capturing the orientational ordering of the entities 
in nematic and/or smectic phases while being a simple 
analytical expression. Hence, it has become one of the 
first choices to model liquid-crystals |3MlO|. 

In the introduction of the BP work it is said: The po- 
tential must have two characteristics: It must be mathe- 
matically simple, involving only functions which are easy 
to calculate; and it must not violate too strongly our sense 
of what is physically correct. Almost forty years later, 
with the huge advance of computing power, the develop- 
ment of very efficient and easy to use software based on 
the so called force-fields, and the possibility of simulating 
a complete and solvated protein by handling several thou- 
sands of atoms, one is tempted to think that these lines 
are obsolete. Notwithstanding, there is nowadays a great 
effort on developing and studying pair po tentials which 
goes in the coarse-grained direction [lll-[l5j . In fact, ellip- 
soids are frequently used as models for highly anisotropic 
and relatively rigid supramolecular structures, such as 
biological membranes, DNA-phospholipidic complexes, 
and other anisotropic entities like nanotubes, inorganic 
nanorods, clay crystals, among a wide variety of colloidal 
particles. Additionally, and hardly contrasting with the 
huge advance of computing power, even for the most sim- 
ple of the pair potentials, i. e. for hard spheres, there are 
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questions that still cannot be answered [laLLZl- Thus, it 
seems that these BP lines are not only current these days, 
but also they will remain valid for many more decades. 

An ellipsoidal shape can be handled to match the shape 
of many molecules and colloids. For this purpose, how- 
ever, it is desirable to know the exact shape of the model 
hard core, as well as its volume and surface. Unfortu- 
nately, the BP hard potential does not provide a defined 
shape and volume. On the other hand, the exact el- 
lipsoidal hard core interaction is only numerically solv- 
able [15|, LL8|, LL9|. Thus, the development of better an- 
alytical approaches to the exact solution is convenient. 
In this regard, Rickayzen [20] proposed a modification 
of the original BP analytical expression (MBP) aiming 
to correct the well-known T-shape mismatch of BP. In 
this work, we slightly modify the Rickayzen expression 
(adding a fitting parameter) to approach even further 
the exact solution. The resulting expression presents a 
good balance between precision and complexity, and so, 
it reasonably fulfills the lines written by BP [1]. Then, we 
tested the analytical expression against the exact solution 
by comparing the corresponding equations of states for 
particles with 1:5 aspect-ratio. For this purpose, we im- 
plemented the replica exchange Monte Carlo simulation 
method. We observed a very good agreement between 
the exact and the MBP analytical equations of state. 
The MBP approach allowed us to study size effects on 
the results reported in ref. [2l| Thus, for the thermody- 
namic limit, we estimate the occurrence of the exact 1:5 
hard ellipsoid isotropic-nematic transition at the volume 
fraction 0.343 ±0.003, and the fluid-solid transition in the 
volume fraction interval (0.592 ±0.006) - (0.634 ±0.008). 
Finally, another crystal structure was captured in coex- 
istence with that found in previous work [21] . 

The paper is organized as follows. The hard ellipsoidal 
pair potential models are given in the next section. There 
we compare the BP and the modified BP predictions with 
those provided by the exact numerical solution. Section 
III describes the replica exchange method for hard bod- 



ies. The equation of state and the structure of the stud- 
ied systems are given in Section IV. Finally, we tackle 
the conclusions in Section V. 



II. HARD ELLIPSOIDAL MODELS 

In this section we summarizes two useful tools to deal 
with hard ellipsoidal interactions. These are the Berne 
and Pechukas analytical approximation and an algorithm 
to produce the exact solution for a given precision. The 
following two subsections present these methods. A third 
subsection is devoted to introduce a modification of the 
Berne and Pechukas analytical expression which improves 
its overall performance. The improvement is shown in a 
final subsection. 

I 



A. Berne and Pechukas 

The Berne and Pechukas potential (BP) [1] is a reason- 
able approximation for the exact interaction of two equal 
hard ellipsoidal particles. It is analytical, mathematically 
simple, easy to implement, fast to compute, and it can 
be used to study the condensed phase of a collection of 
prolate or oblate particles via numerical experiments |22J - 
[25j. In this approach, molecules are represented with an 
uniaxial ellipsoidal Gaussian and their interaction is then 
related to their overlap integral. In this way, a coarse- 
grained potential is built which successfully captures the 
anisotropic nature of the entities. The expression for the 
distance between the geometric centers of the ellipsoids 
when the particles are at contact, cr#p(ui, iij, r), is given 
by 



a B p(ui,nj,f) = <j_l( 1- -x 



(r • u^ + r • Uj) 2 (r • u^ — r • Uj) 2 



XUi • u^ 



X Ui • Uj 



-1/2 



(1) 



where u^ and Uj are the versors (unit vectors) along the 
main axis of each particle, and r is the versor along the 
line joining the geometric centers. Here, the anisotropy 
parameter \ is 



X : 



(2) 



where cry and a± are the parallel and perpendicular di- 
ameters with respect to the main axis, respectively. 



B. Exact numerical solution 

An ellipsoidal surface centered at r^ and oriented ac- 
cording to iii (parallel to its main axis) is given by the 
following quadratic form 



Ai(r e ) = ( r e - n Y • Ai • ( r e - r» ), 



(3) 



where the scalar Ai(r) = const, and r e is a point at the 
surface. In particular, the ellipsoid surface is represented 
by Ai(r e ) = 1, having 



A i= U*(ui) 



■U(ui), 



(4) 



the rotation matrix U(u^) (the matrix which converts a 
vector from the space-fixed to the body-fixed coordinate 
system |26j), its transpose U*, and the diagonal matrix 



(5) 




r 



Thus, the geometry of the particle is given by D and its 
orientation is given by U(ui) (or simply by Uj). 

Let's consider two equal, arbitrarily oriented, and non- 
overlapping ellipsoids i and j, at contact at point r c . The 
vector normal to the i surface at the contact point is 



ni(r c ) = VA(r c ) 



( 



)• 



(6) 



A similar equation can be written for the vector normal 
to the j surface Uj. Since the tangent plane is common 
for both ellipsoids, the normal versors n^ = n^ / 1 n^ | and 

n J = n j/\ n j\ fulfill 



n*(r c ) + n j (r c ) = 0. 



(7) 



This condition was originally employed by Perram and 
Wertheim for developing an algorithm for numerically 
determining the point r c yj| [l9| . In their work the El- 
liptic Contact Function (ECF) is introduced, which con- 
tains the information given by equation (0 and allows 
determining the distance of closest approach. Later, the 
ECF procedure is reviewed by Paramonov and Yaliraki, 
who contributed with a clear geometric interpretation of 
the Perram and Wertheim approach [15]. The expres- 
sion for the function that connects the particles centers 
through the geometric place where the vectors V^4i(x c ) 
and VA/(x c ) are antiparallel is given by [15| 



= (A) 



+ (l-A)A i ) 1 -(AA i -r i + (l-A)A j -r j ) (8) 



where A G [0, 1] is a scalar parameter. Note that for A = 1 
and the geometric centers of the ellipsoids i and j are 
obtained, respectively. The contact point r c lies on this 
trajectory and corresponds to a unique value of A, A c , 
which fulfills A c G (0, 1). Furthermore, *4*(r c ) = A/(r c ), 
with r c = x c (A c ). 



With the above expressions it is easy to implement an 
iterative procedure to yield r c with the desired precision. 
In particular, we implemented the bisection algorithm. 
It starts by evaluating A (A) = ^(x c (A)) — A/(x c (A)) 
for A = 0.5 (A(A) is a monotonously decreasing function 
and is zero solely at A c ). A positive A (A) means 1 > 
A c > A and so, A is increased in such a way to reduce 
in half the interval. Conversely, a negative A (A) means 
< A c < A and A is decreased accordingly. In this way 
the interval is reduced as l/2 n , being n the number of 
iterations. This method is simple and safe but not the 
faster. Approximately 20 iterations yields an error of A 
smaller that 1 x 10 -6 . Note that the involved operations 
are only products and summations which translate into 
a fast computation. 

Additionally, the contact parameter A c defines the ex- 
treme value of the linear combination of the quadratic 
forms Ai(~x c ) and *4j(x c ), namely 

5(A) =XAi (x c (A)) + (1 - X)A 5 (x c (A)) . (9) 

In other words, < 5(A) < <S(A C ) (<S(A) is an strictly 
concave function which ensures it has a unique maxi- 
mum at A c ). Note that from equations (J3j), (j4]), and (J9j), 
5(0) =5(1) =0. The contact parameter also defines the 
Perram-Wertheim (PR) contact distance <Jpw, 



and 



&PW 



V^c 



(10) 



Hence, ellipsoids having their geometric centers sepa- 
rated a distance r smaller than cjpw overlap whereas 
they do not for r > crpw- In other words, (5(A C )) -1 / 2 
can be interpreted as a scaling factor needed to bring the 
particles into contact. In particular, the BP analytical 
expression for the contact distance corresponds to [15] 



°BP 



VW) 



(ii) 



which, in general, overestimates the exact result, app > 
crpw, since 5(A C ) > 5(1/2). 



C. Modified Berne and Pechukas 

The BP approximation to the contact distance of two 
equal ellipsoids becomes poor for Uiiij — )> 0. In particu- 
lar, for case E of Fig.CQthe BP result is ((a^ +a|)/2) 1 / 2 
instead of (a± + crii)/2, which is the exact solution. To 
improve the overall behavior of the BP approximation 
Rickayzen proposed [20] 



*mbp = <T±{l-\x[A + +A-] + {l-x)x'[A + A-]~ f 



being 



(f • Uj ± f • Uj 
1 ± XUi ■ iij 



-1/2 

(12) 
(13) 



X 



cr± 



cr± 



(14) 



In his work 7 is unity. These expressions lead to very 
good results, since the introduced term does not alter 
the correct answers given by BP for cases A-D and F of 
Fig. [1] (note that in these cases A+A~ = 0). Additionally, 
x' is given such that case E is also satisfied. Thus, the 
above expressions correctly describe all cases shown in 
Fig. [TJ Also note that (Jmbp < &bp for all aspect ratios 
and 7. 

There are, however, other expressions which produce 
the correct result of cases A-F of Fig. [TJ After trying 
some of them, we concluded that expressions ([T2]) - (fT4]) 
(MBP) show a good compromise between accuracy and 
complexity. We solely added 7 as a free parameter to be 
adjusted as a function of the aspect ratio. 



D. Comparing BP and MBP to the exact 
numerical solution 



In the previous section a modification of the BP (MBP) 
is presented, which fixes the T-shape mismatch while pre- 
serving the goodness of the BP analytical solution for the 
rest of the cases. We now compare the BP and the MBP 
predictions with the corresponding exact numerical solu- 
tion (the error of the numerical solution is always below 
10 -8 ). This is shown in detail for the 1:5 oblate case in 
Fig. [2] a). This plot shows the appearance frequency of 
a given deviation percentage as a function of the devi- 
ation percentage, e» = 100(a arii - cr PWi )/a PWi . Here, 
a ani is the analytical prediction of the overlap distance 
for the random position and orientation z, and crpWi is 
the corresponding exact solution. The plot is then built 
by considering 7V C = 1 x 10 8 different positions and orien- 
tations (note that for a spatially fixed ellipsoid we must 
generate a random position, keeping r fixed, and a ran- 
dom orientation of the second ellipsoid to sample from 
the space of all possible configurations). The black, cyan, 
and dashed red lines represented in the plot correspond 
to the BP and the MBP results with 7 = 1 and 7 m i n , re- 
spectively. Here, jrnin is the value of 7 which minimizes 
E — ^/N c ^2 i c \ei\. As can be seen, the BP deviations 
are greater or equal than zero, meaning that BP overes- 
timates the volume of the corresponding oblate ellipsoid. 
Additionally, a large peak at zero is seen, which indicates 
that many configurations are correctly described. In the 
same direction, we should mention that 50% of the ran- 
domly generated configurations have a deviation smaller 
than 3.71%. Nonetheless, the distribution also shows a 
long tail reaching deviations as large as 20%. A look to 
the corresponding snapshots for large deviations shows 
that the ellipsoids are close to the T-shape configuration 
(not shown). 

Fig. [2] a) shows the main goal of the MBP, i. e., the 
long tail towards positive deviations of the BP expression 
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FIG. 1. Two ellipsoids at different contact configurations. The exact center-center distance at contact, a, is given as a 
function of cry and cr±. The relations among the normal versors u^, Uj and the center-center vector r are also shown. The BP 
approximation produces the exact solution for all cases but E. MBP yields the exact solution for all shown cases. 




FIG. 2. a) Appearance frequency, F, against the percentage 
of deviation of the analytical overlap distance with the cor- 
responding exact distance, e% — 100(cr a n i — crpw i )/o'pw i , for 
the 1:5 oblate case. The histogram is built by considering 
N c = 1 x 10 8 random positions and orientations, b) Average 
deviation, E = ^-/N c ^ i c \a\ (bullets), and the optimized 
7 parameter, j m in (open symbols), for several aspect ratios. 
The vertical dark line at 1:5 points out the case shown in 
panel a). Black solid lines, cyan solid lines, and red dashed 
lines correspond to the BP, MBP with 7=1, and MBP with 
jmin analytical expressions, respectively. 



is strongly suppressed. This occurs for both, 7 = 1 and 
7 = 7mm- That is, for all considered cases, we did not de- 
tect deviations larger than 6.04% and 4.04%, respectively. 
From the simulation point of view this is a nice result 
since it guaranties configurations with r < 0.959ctmbp 



to be overlapped. In the same line of thinking, config- 
urations having r > gbp cannot be overlapped. These 
results aid reducing the number of configurations to nu- 
merically solve when the purpose of the simulation is to 
work with the exact hard core model of the ellipsoids. 

Fig. [2] a) also shows that a price must be paid when 
correcting the BP T-shape mismatch. That is, MBP in- 
troduces negative deviations. Consequently, there are 
overlapped configurations that MBP will consider as free 
of overlaps. However, we did not detect these deviations 
to exceed 0.84% and 1.39% for 7 — 1 and 7 = 7 m i n , re- 
spectively. Thus, ellipsoids are definitely not overlapped 
for configurations having r > 1.020ctmbp- On the other 
hand, the average deviation, E — ^/N c ^ ji c |e^|, is re- 
duced from 5.41 to 1.13 (7 = 1) and to 0.62 (7 = j m in) 
when implementing MBP instead of BP, which points out 
that the true ellipsoidal shape is much better approached. 

In brief, Fig. [2] a) shows an improvement of MBP 
when compared to BP for the 1:5 aspect-ratio oblate 
case. Now, what happens when considering other aspect 
ratios or prolate ellipsoids? The answer to this ques- 
tion is given in Fig. [2]b). There, the average deviation, 
E = l/iV c ^ c |ei|, as a function of the aspect ratio is 
shown for both, oblate and prolate ellipsoids. Again, the 
black line is used to represent the BP results whereas 
the cyan line and the red dashed line correspond to the 
MBP expression with 7 = 1 and 7 = 7 m i n , respectively. 
The vertical dark line at 1:5 points out the case shown 
in Fig. [2] a), i. e., the 1:5 aspect-ratio oblate case. This 
panel shows that the MBP expression approaches bet- 
ter the exact case for all studied aspect-ratios -from the 



1:25 oblate (leftmost points) to the 25:1 prolate (right- 
most points)- and for both, 7 = 1 and 7 = 7 m i n . The 
only exception is the 1:1 ellipsoid (sphere) where all ap- 
proaches yield the exact result. Prolate ellipsoids, how- 
ever, are reasonably described by BP for all aspect-ratios, 
and thus, the gain by implementing MBP in these cases 
is not very large. Conversely, the main BP deviations oc- 
cur for oblate shapes reaching E values over 10% for the 
1:20 and 1:25 cases. In these cases the gain of MBP is 
remarkable, since it always shows E values below 4%. In 
view of these results, it seems that the gain in precision 
is worth the little extra operations needed to compute 
MBP instead of BP. In addition, the gain of employing 
7 = 7min instead of 7 = 1 is close to 50% (see Fig. [2] 
b)). The values of j m in are shown in Fig. [2]b) as open 
symbols. For 1:5 oblates j m i n = 0.794. 



III. REPLICA EXCHANGE MONTE CARLO 

Even though the analytical expressions for determin- 
ing whether or not two ellipsoids overlap are relatively 
fast to compute, sampling from crowded systems is al- 
ways a difficult task [17]. We expect 1:5 oblate hard 
ellipsoids to show a fluid-nematic transition at relatively 
low densities [27-29] and a nematic-crystal transition at 
high densities [27|, [28[ , and this work attempts to cap- 
ture both. By analogy with the hard sphere case, we also 
expect the crystal equilibrium branch of the phase dia- 
gram to lie below a metastable nematic phase branch (the 
continuation of the nematic branch toward larger densi- 
ties than the crystallization point), which surely hinders 
equilibrating the system at high densities [30| . For these 
reasons, we are implementing the replica exchange Monte 
Carlo methodology, which is well proven to assist the sys- 
tems to reach equilibrium at difficult (high density / low 
temperature) conditions. 

In the classical replica exchange scheme, n r identical 
replicas are considered each following a typical canoni- 
cal simulation at different temperatures |3ll-l33|. Thus, 
an extended ensemble can be defined so that its parti- 
tion function is Q extended = O^i Qnvti, being Qnvti 
the partition function of ensemble i at temperature T^, 
number of ellipsoids TV, and volume V. The existence 
of this extended ensemble justifies the introduction of 
swap trial moves between any two ensembles (each be- 
ing sampled by only one replica at a time), whenever the 
detail balance condition is satisfied. If all (i,Ti)(j,Tj)—)> 
(j,Ti)(i,Tj) swap trials have the same a priori probabil- 
ity of being performed, the swap acceptance probability 
becomes 

P acc = min(l,exp[(P j -p i )(U i - U 3 )}) (15) 

where ft = l/^ksTi) is the reciprocal temperature of 
replica z, ks is the Boltzmann's constant, and U{ is the 
energy of replica i. Hence, by introducing these swap 
trials, a particular replica travels through many temper- 
atures allowing it to overcome free-energy barriers. Ad- 



ditionally, sampling on particular ensembles is not dis- 
turbed but enriched by the different contributions of the 
n r replicas. 

Since we are dealing with hard ellipsoids the tempera- 
ture plays a trivial role. Then, to take advantage of the 
method, we must perform the expansion in pressure in- 
stead of temperature [30]. Hence, the partition function 
of the extended ensemble is given by [30, l34j 



/extended 



=n< 



{NT Pi 



(16) 



where QNTPi is the partition function of the isobaric- 
isothermal ensemble of the system at pressure P^, tem- 
perature T, and with TV particles. 

This extended ensemble is sampled by combining stan- 
dard NTPi simulations on each replica (involving trial 
displacements, rotations of single ellipsoids, and trial vol- 
ume changes) and replica exchanges (swap moves at the 
replica level). To satisfy detailed balance, these swap 
moves are performed by setting equal all a priory prob- 
abilities for choosing adjacent pairs of replicas and using 
the following acceptance probability [30] 

P acc = min(l,exp[/3(P i - Pj)(Vi - V s )]), (17) 

where V{ — Vj is the volume difference between replicas 
i and j. Adjacent pressures should be close enough to 
provide reasonable exchange acceptance rates between 
neighboring ensembles. In order to take good advantage 
of the method, the ensemble at the smaller pressure must 
also ensure large jumps in configuration space, so that the 
larger pressure ensembles can be efficiently sampled. 

The probability for selecting an ellipsoid displacement 
trial, Pd, an ellipsoid rotation, P r , for selecting a volume 
change trial, P V: and a swap trial, P s , are fixed to 



Pd 
P v 

Ps 



P r =n r N/(n r (2N + l) 
n r /(n r (2N + l)+w), 
w/(n r (2N + l)+w), 



•w), 



(18) 



where w <C 1 is a weight factor. Note that Pd + P r + Pv + 
P s = 1, as it should. The probability density function to 
have the next swap trial move at the trial n t is given by 



P{n t ) = P s exp(-P s n t ). 



(19) 



Hence, we may obtain the next swap trial move from 
n t = — ln(£)/P s , with £ being a random number uni- 
formly distributed in the interval (0,1) [35|, |36[. We set 
all ellipsoids of a given replica to have the same a priori 
probability of being selected to perform a displacement 
or a rotational trial. The same is true for selecting a 
replica to perform a volume change trial. 

The trials [l,n t — 1] are displacements, rotations, and 
volume changes, and so, they can be independently per- 
formed on the replicas. This has the advantage of being 
easily parallelized. The algorithm is parallelized through 
message passing interface (mpi) fortran in n r threads, 



though quad core desktops are used. Since all swap tri- 
als are performed in a single thread, the efficiency of the 
parallelization increases with decreasing w. We employed 
w = 1/100. Verlet neighbor lists [37|,l38| are used for sav- 
ing CPU time, which can be quite large for the replicas 
evolving with the highest pressure values. These lists 
must take into account displacements and rotations to 
timely update. 

Simulations are performed in two steps. All simula- 
tions are started by randomly placing the ellipsoids in a 
random orientation (avoiding overlaps), so that the ini- 
tial volume fraction is cp = v e p = 0.2, where p is the 
number density, v e = 47r<7||<7 2 L /3 is the ellipsoid volume 
(for all models), an =1, and a_\_ = 5. We first perform 
about 2 x 10 13 trial moves at the desired state points, dur- 
ing which we observe that the replicas reach a stationary 
state (equilibrating procedure). We then perform 2 x 10 13 
additional trials during which various measurements are 
carried out, with results described in the following sec- 
tion. 

The maximum particle displacements, maximum rota- 
tional displacements, and volume changes for trial moves 
are adapted for each pressure to yield acceptance rates 
close to 0.3. Thus, particle displacements, rotations, and 
volume changes of ensembles having high pressures are 
smaller than those associated to ensembles having low 
pressures. An optimal allocation of the replicas should 
lead to a constant swap acceptance rate for all pairs of 
adjacent ensembles. For a temperature expansion, the 
efficiency of the method peaks at swap acceptance rates 
close to 20% [39| . In this regard we implemented a simple 
algorithm to smoothly adjust the state points (pressures) 
to yield an approximately constant swap acceptance rate 
while keeping constant the maximum and minimum pres- 
sure. At the starting point, we use a geometric progres- 
sion of the pressure with the replica index. Note that 
the adaptation of maximum displacements and the shift 
of the desired pressures violates the detail balance con- 
dition. Thus, these procedures are performed only dur- 
ing the equilibrating procedure (the first 2 x 10 13 trial 
moves). This work is performed by considering N = 100 
and N = 200 ellipsoids, and n r = 64 (to cover a wide 
range of densities while keeping swap acceptance rates 
over 20%). 



IV. RESULTS 

Results are split in three subsections. The first one 
is devoted to compare the equation of state obtained 
by considering the exact model of the 1:5 oblate el- 
lipsoids with the analytical expression of Berne and 
Pechukas (BP) and its modification (MBP) presented in 
subsection III CI This is done by considering 100 particles 
(N=100). The second subsection compares the N=100 
MBP case with 7 = jmin to simulations performed with 
N=200. The obtained structures are analyzed in a third 
subsection for this larger system size. 




FIG. 3. a) Equations of state of the exact and analytical hard 
1:5 oblate ellipsoidal models, Z(ip) (if is the volume fraction), 
b) Isothermal compressibility obtained from density fluctu- 
ations, x(^)- c ) Order parameter, Qe(cp). For all panels, 
black circles, cyan squares, red crosses, and blue plus symbols 
correspond to the exact PW overlap distance, and the analyt- 
ical solutions given by BP, MBP with 7=1, and MBP with 
7 = 7min, respectively. All data correspond to N = 100. 



A. Exact vs Analytical overlap distance 

The exact hard 1:5 oblate ellipsoidal model can be 
studied for a moderate system size. This is possible since, 
on the one hand, the exact iterative procedure to solve 
the apw overlap distance is relatively fast, and on the 
other, the numerical solution is computed only for those 
few cases where the ellipsoids geometric centers distance, 
r, is less than &bp- In this way the exact PW simu- 
lation is approximately three times more computation- 
ally demanding than the analytical cases. As mentioned, 
n r = 64 replicas are employed to cover a wide pressure 
range while keeping large swap acceptance rates. The 
probability distribution functions (PDFs) obtained for 
this system are given in previous work [21]. From them 
the dimensionless pressure Z = (iP/p and the isother- 
mal compressibility \ are calculated as a function of the 
most frequent volume fraction, ip. These functions are 
given in panels a) and b) of Fig. [3] as black circles. The \ 
values are obtained by means of the density fluctuations, 
i. e., by x = N(< p 2 > - < p > 2 )/ < p > 2 , which 
should equal \ — dp/d(PP) according to the fluctuation- 
dissipation theorem. Finally, panel c) shows the order 

parameter Q 6 = (ff YZZ-e I <*6m(0, <t>) > l 2 )^ be " 
ing < YsmiO^cj)) > the average over all bonds and con- 



figurations of the spherical harmonics of the orientation 
polar angles and <fi |30l 140, l4lj . Qe approaches zero for 
a completely random system of a large number of points, 
and increases when configurations present angular order. 
All these data shown as black circles correspond to the 
exact solution and are taken from ref. [2l| . Also from this 
reference, we are including the data corresponding to the 
BP overlap distance approximation as cyan squares. 

Simulations for the exact PW case show a isotropic- 
nematic transition at cp = 0.341 and a nematic-crystal 
transition at the volume fraction interval of 0.584 — 
0.605 [21]. The fluid-fluid transition is evidenced by a 
plateau of Z(tp), a kink of x(<p), and a practically invari- 
ant Qq((p). The fluid-crystal transition is characterized 
by a jump of ip accompanied with a kink of x(^) and a 
steep increase of Qs(ip). For the isotropic region Z(ip) 
agrees with the simulation data provided by Mc. Bride 
and Lomba [42[ and so, our Z(cp) is describe d by the Vega 
equation of state for the isotropic fluid [42|, |43[ . The BP 
model also captures the three phases, although the tran- 
sitions are shifted to lower densities, and the low and high 
density branches differ from the exact case (see the insets 
of Fig. [3]). The transitions occur at cp = 0.274 and at the 
range 0.577 — 0.595. These discrepancies are due to the 
fact that gbp > crpw For the nematic phase, however, a 
very good agreement between BP and PW is seen for all 
panels of Fig.[3j This is due to gbp well approaches crpw 
for parallel configurations (see Fig. []]). Discrepancies are 
more pronounced at low densities where T-shape configu- 
rations are frequent, for which gbp may reach 1.2 x crpw 
(see Fig. [2]). Thus, the MBP approach, which corrects 
the T-shape BP mismatch, is expected to approach bet- 
ter the PW simulations. 

The MBP data with 7 = 1 are shown as red crosses 
in the three panels and the insets of Fig. [3l As can be 
seen, discrepancies between the exact and the analytical 
cases are clearly diminished. In particular, the isotropic- 
nematic transition is now observed at cp = 0.331, and 
so, the relative difference with the exact case diminishes 
from 20% to 3%. Notwithstanding, the exact-analytical 
agreement can be further improved by setting 7 = jrnin- 
The corresponding results are also included in Fig. [3] as 
blue plus symbols. In this case all discrepancies practi- 
cally vanish (though a somewhat more structured crystal 
is still produced according to Qq). The obtained transi- 
tions occur at tp = 0.344, and in the range 0.581 — 0.599. 



B. Size effects 

The faster computation of the analytical case allows 
us to study a larger system size, i. e., N = 200. For 
this system size and keeping n r = 64 we obtain swap 
acceptance rates close to 30%. These high acceptance 
rates are a consequence of the large overlapping areas 
of the probability distribution functions (PDFs) shown 
in panel a) of Fig. |4j There, the overall trend of the 
PDFs to get narrower and higher with increasing den- 



^0.03 

Q 

Ph0.02 




FIG. 4. a) Probability distribution functions (PDFs) of vol- 
ume fraction fluctuations for each of the n r — 64 pressure 
values and for N = 200. b) Equations of state, Z(ip). c) 
Isothermal compressibilities, x(^)) obtained from density fluc- 
tuations, d) Order parameters, Qe((p). For panels b), c), and 
d), dark circles and light squares correspond to N = 200 and 
N — 100, respectively. Vertical dotted lines highlight the 
N — 200 transitions. All data correspond to the analytical 
MBP case with 7 = ~fmin ■ 



sity (pressure) is seen. This makes a general decrease 
of the isothermal compressibility x — dp/d(PP) evident. 
The overall trend is disrupted at the transitions where 
PDFs are distorted from their natural Gaussian shape. 
Indeed, PDFs turn bimodal at the nematic-crystal tran- 
sition. In brief, all features captured by the N = 100 sys- 
tem also appear for the N = 200 case (panel a) should be 
compared to Fig. 1 of previous work [2l|). Nonetheless, 
some differences appear: The PDFs are higher and nar- 
rower due to the larger system size, the nematic-crystal 
transition density gap enlarges and shifts to larger den- 
sities, and a small disruption of the PDFs heights trend 
appears at very large ip (ip ~ 0.65). In particular, the 
nematic-crystal transition shift is clearly shown in panels 
b), c), and e) of the same figure, where the N = 100 
data are included as cyan squares for an easy compari- 
son. This is the only practical difference of Z(ip). x(^) 
evidences that size effects are absent from the isotropic- 
nematic transition (this is in agreement with Allen and 
Mason results for the 3:1 prolate case kj and opposed to 
the Zarragoicoechea et al. suggestion |45j). On the other 
hand, it reflects the small disruption of the PDFs heights 
trend as a little kink at high densities. Finally, Qq{^p) 
suggests the formation of less defective crystals for the 



larger system size (Qq is larger for N = 200 at densities 
above the nematic-crystal transition). 

For N = 200, the isotropic- nematic transition occurs 
at (f = 0.345 and Z = 8.6. The nematic-solid transition 
is located in the volume fraction interval 0.585 — 0.614 
for Z = 27.7. Thus, we can now extrapolate the val- 
ues corresponding to N = 100 and 200 to N — >> oo. 
This procedure leads to cp = 0.346 and Z = 8.7 for 
the isotropic-nematic transition and to ip in the range 
of 0.589 - 0.629 at Z = 29 A for the nematic-solid tran- 
sition, which are our estimates for the MBP model at 
the thermodynamic limit. Since a similar shift is ex- 
pected for the exact case, our estimate for the ther- 
modynamic limit for the isotropic-nematic transition is 
if = 0.343 ± 0.003 and Z = 8.4 ± 0.3, and in the vol- 
ume fraction interval (0.592 ± 0.006) - (0.634 ± 0.008) at 
Z = 30.4 ± 0.9 for the nematic-solid transition. These 
values can be compared to those shown by the phase di- 
agrams given in refs. (27|, [28j. We estimated them to be 
ip ~ 0.37 for the isotropic-nematic transition and the 
interval (p ~ 0.60 — 0.67 for the nematic-solid transi- 
tion. Additionally, Samborski et. al. placed the isotropic- 
nematic transition at the ip range of 0.333 — 0.351 [29]. 
Hence, our value (taken where function x{v) peaks) lies 
exactly in the middle of this range. We should also 
add that replicas with ip ~ 0.343 are not totally ne- 
matic or isotropic according to the observed snapshots 
(not shown). Thus, the transition takes place over a vol- 
ume fraction range instead of a single point. Taking the 
two intersection points of x( ( / ? ) with the horizontal line 
that contains the local minimum of x(^) ( a ^ P ~ 0.312 
for N = 100 and the exact model) we obtain the vol- 
ume fraction range 0.312 — 0.363. This range translates 
into 0.314 — 0.365 by applying the corresponding shift to 
estimate the transition at the thermodynamic limit. 



C. Structure 

Fig. [5] shows the radial distribution functions, g{r), 
and order parameter, p(r) =< l/2(3(u^ • Uj) 2 — 1) >, for 
different pressures. Panel a) corresponds to the lowest 
pressure where both radial functions signal an isotropic 
phase. That is, the g(r) smoothly increases peaking at 
r ~ <jx, while the p(r) shows no sign of long range an- 
gular order at large distances (the p(r) peak at r = cry 
results from the few and forced parallel configurations 
at short distances). The corresponding snapshot clearly 
shows this fact. Panel b) corresponds to the structure 
found for cp ~ 0.45. In this case the function p(r) shows 
a long range alignment of the ellipsoids though the g(r) 
still points to a fluid phase (nematic). The only peak 
shown by the g(r), as in the previous case, corresponds 
to r ~ a±, but now the peak is higher probably due 
to the side-to-side configurations in the developed lay- 
ers (see the corresponding inset). The last two panels 
correspond to a solid phase. Pressure in c) is close to 
the transition and pressure in d) is the highest. Both 




FIG. 5. Radial distribution functions (black lines), g(r), and 
their corresponding order parameter (cyan lines), p(r), for 
different pressures (increasing from left to right). Panels a), 
b), c), and d) correspond to the isotropic, nematic, crystal A, 
and crystal B structures. The insets show the corresponding 
snapshots. 



crystals differ from each other. The structure of panel 
c) (crystal A) is the one found for the TV = 100 case 
of previous work [21]. It is characterized by a first and 
large peak at r ~ 1.2, which corresponds to the touch- 
ing and stacked ellipsoids shown in the inset. There is a 
second peak at r ~ 2.4, corresponding to two ellipsoids 
separated by a third one and sandwiched by them (all 
belonging to the same stack). Finally, the wide peak at 
r ~ 4.3 corresponds to side-to-side configurations of par- 
ticles belonging to different stacks. This peak is smaller 
than a± since the particles of the adjacent stacks are par- 
tially sandwiched. The structure of panel d) (crystal B) 
shows a shorter main peak which is also shifted to r ~ 1.5 
(the shoulder at r ~ 1.2 and the small peak at r ~ 2.4 
correspond to contributions to the average from replicas 
having a crystal A like structure). In this arrangement, 
the secondary peak appearing at r ~ 2.4 is also shifted to 
r ~ 3.0. Thus, in both cases, the second peak appears at 
a distance two times larger than that of the main peak. 
Other differences appear at larger distances which are 
related to the way the stacks are arranged. 

The differences found for the radial distribution func- 
tions can be understood from the structural details shown 
in figure El From this figure it is seen that the projection 
of the particles geometric centers towards a plane per- 
pendicular to the stacks axes are arranged in a hexago- 
nal lattice for crystal A (left panel), whereas they pro- 
duce a square centered pattern for crystal B (right panel) . 
Additionally, the intra-stack arrangements also differ re- 
spect to the other. In crystal A the entities of a given 
stack are more perpendicular to the stack axis than for 
the crystal B structure. In other words, the angles 
a =< 90° - cos _1 (u* • r 8 ) > are 62° and 44° for crystal A 
and crystal B, respectively. The more tilted intra-stack 
arrangement of crystal B leads to an increase of the num- 
ber of stacks to be arranged for a given area (perpendic- 
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FIG. 6. Detail of the different crystal structures found in the 
N — 200 system (cases d) and e) of Fig. 2J. The shown a 
angles, defined as < 90° — cos -1 (iii • r s ) > where r s is the 
versor pointing along the stack axis, are close to 62° and 44° 
for crystal A and B, respectively. 



ular to the stacks' axes), though the stacks contain less 
particles for a given axis depth. This more tilted intra- 
stack arrangement explains the shift of the first two peaks 
of the g(r) towards larger distances of crystal B. 

Note that crystal B is not detected for a system size of 
N = 100. Additionally, the crystal A structure found for 
N = 100 shows many more defects than the one found 
for N = 200. This is evidenced by the smaller Qq values 
and also directly seen from the snapshots. Hence, size 
effects are important for determining the correct crys- 
tal phase (implementing non-orthogonal unit cells by al- 
lowing the lattice vectors to change may improve this 
point). Accordingly, we can only state that the structures 
here reported correspond to the studied system size. For 
N = 200, we observe that replicas having a structure 
like crystal B are preferably located at higher pressures, 
though more replicas produce the crystal A structure. 
Consequently, a kink appears in xiv) a ^ high pressures, 
separating a rich crystal A region at smaller pressures 
from a rich crystal B region at higher pressures. This 
points out that crystal B has a larger compressibility 
than crystal A and may suggest a crystal-crystal tran- 
sition. Finally, we should also mention that the struc- 



tures here reported are different from those pointed out 
by Donev et. al. [46|, |47j. They reported that structures 
composed by a two-layer lamination where layers are dis- 
posed rotated ir/2 with respect to the other yields a pack- 
ing fraction of (p = 0.7707. To obtain this large packing 
fraction layers must be face-centered square planar and 
have all ellipsoids oriented with one of their semiaxes 
perpendicular to the layer and the other two oriented 
along the axes of the lattice. The authors also mention 
that there is nothing suggesting that this family of struc- 
tures is the densest. On the other hand, a fit of the 
form Z ~ (jp d — Lp)~ x to the high pressure branch of 
Z((p) for TV = 200 (panel b) of Fig. gj) leads to a diver- 
gence at (fd = 0.769 ± 0.002, which is virtually equal to 
(f = 0.7707. Consequently, several crystal candidates for 
the equilibrium structure at high densities seem to exist. 



V. CONCLUSIONS 

We used the replica exchange Monte Carlo technique 
to produce the equation of state of hard 1:5 aspect-ratio 
oblate ellipsoids for a wide density range. In addition 
to the analytical approximation of the overlap distance 
given by Berne and Pechukas (BP) and the exact numer- 
ical solution given by Perram and Wertheim, we imple- 
mented a simple modification of the original BP approx- 
imation, which corrects the known T-shape mismatch of 
BP. We found that this approximation produces an equa- 
tion of state practically equal to that obtained for the 
exact overlap distance solution. We then used this ap- 
proach to study a larger system size (N = 200). The 
produced results allowed us to estimate the locations of 
the isotropic-nematic and nematic-crystal transitions for 
the thermodynamic limit at ip = 0.343 ± 0.003 and Z = 
8.4±0.3 and in the interval (0.592±0.006)-(0.634±0.008) 
for Z = 30.4 ± 0.9, respectively. 
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